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Abstract 



We present an analytical continuous equation for the Tang and Leschhorn 
model [Phys. Rev A 45, R8309 (1992)] derived from his microscopic rules 
using a regularization procedure. As well in this approach the nonlinear term 
(V/i) 2 arises naturally from the microscopic dynamics even if the continuous 
equation is not the Kardar-Parisi-Zhang equation [Phys. Rev. Lett. 56, 889 
(1986)] with quenched noise (QKPZ). Our equation looks like a QKPZ but 
with multiplicative quenched and thermal noise. The numerical integration of 
our equation reproduce the scaling exponents of the roughness of this directed 
percolation depinning model. 

PACS numbers: 47.55. Mh, 68.35.Fx 
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The investigation of rough surfaces and interfaces has attracted much attention for 
decades due to is importance in many fields, such as the motion of liquids in porous media, 
growth of bacterial colonies, crystal growth, etc. When a fluid wet a porous medium, a 
nonequilibrium self-affine rough interface is generated. The interface has been characterized 
through scaling of the interfacial width w — ([hi — (hi)] 2 ) 1 / 2 with time t and lateral size L. 
The result is the determination of two exponents, f3 and a called dynamical and roughness 
exponents respectively. The interfacial width w ~ L a for t ^> t* and w ~ t 13 for t < f, 
where t* = L a ^ is the crossover time between this two regimes. Much effort has been 
done to understand the leading mechanisms of these processes and to try to explain how 
the dynamics affects the scaling exponents The formation of interfaces is determinated 
by several factors, it is very difficult to discriminate theoretically all of them. The knowl- 
edge of the dynamical nonlinearities, the disorder of the media, and the theoretical model 
representing experimental results are difficult to overcome due the complex nature of the 
growth. The disorder affects the motion of the interface and leads to its roughness. Two 
main kinds of disorder have been proposed: the "annealed" noise that depends only of time 
and the "quenched" disorder due to the inhomogeneity of the media where the moving phase 
is propagating. The discrete models provided an useful approach to obtain the exponents 
that allows its classification in universality classes. By extensively studying these models, 
one can obtain the scaling behaviors and the corresponding universality classes, and then 
associate the continuum stochastic equations with the given discrete growth models. 

The most used method of establishing the correspondence between a continuum growth 
equation and a discrete model, is to numerically simulate the model and compare the ob- 
tained scaling exponents with those of the corresponding continuum equation. In this context 
attempts are being made to classify quenched disorder models in terms of universality classes 
based on equation of motion such as 

= T + uV 2 h + A (v/i)2 + h) + ^ f) (1) 

where T is the driving force responsible of the advance of the interface, £ (x, h) is the 
quenched disorder or pinning forces and ?7(x, t) is the thermal noise. The noises are white. 
Eq. (HD is Kardar-Parisi-Zhang equation [0 with quenched noise (QKPZ). When A — > 
the quenched Edward- Wilkinson equation (QEW) || is recovered. In absence of quenched 
noise his thermal versions are recovered, named KPZ and EW respectively. Much effort 
has been done in order to classify discrete models and experiments in universality classes 
depending on the value of the coefficient A associated with the nonlinearity of the QKPZ. 
Numerical mj studies indicate A is relevant at the depinning transition for discrete models 
in anisotropic media. These results only shows that the nonlinear term exist but they do not 
confirm that these models are represented by the QKPZ. However the exponents obtained 
by numerical simulation of Eq. P), without thermal noise 0, agree very good with those of 
the model in anisotropic media. 

A powerful method of establishing the correspondence between a continuum growth equa- 
tion and a discrete model is to derive the continuum equation from a given discrete model 
analytically. Among them, a systematic method proposed by Vvedensky et al. |7|, where 
the continuum equations can be constructed directly from the growth rules of the discrete 
model based on the master-equation description, has been applied to the derivation of growth 
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equations for some discrete models with thermal noise. This method has proved to be 
useful to derive continuous equation from the microscopic ones with the advantage that the 
sources of the terms of the Langevin equation can be identified and their magnitude related 
to the microscopic dynamics. While the method described by theses authors can be applied 
to general case it is more easy to achieve the same results using a microscopic equation based 
on rules for the evolution of the height. The derivation of continuous equations from discrete 
models is an interesting subject which has not been addressed in the context of growth in 
presence of quenched media. 

The aim of this work is to obtain the continuous equation from the microscopic dynamics 
of a model with quenched noise in order to establish if it is related in some way with the 
QKPZ equation. The main goal of our paper is to have obtained analytically the differential 
equation which describe the dynamics of the Tang and and Leschhorn (TL) model TO]. To 



our knowledge this is the first time that a Langevin equation is obtained from the microscopic 
dynamics in quenched media. 

As we shall show below the dynamics of the height is strongly affected by a multiplicative 
quenched noise. We chose this model because it presents the principal features of some 
experiments like the imbibition of a viscous fluid in a porous media driven by capillary 
forces P"T| , |r2| |. In the TL model, the interface growth takes place in a square lattice of edge 



L with cells of size a that represents the mean size of a pore. For a given applied pressure 
p > , we can divide the cells into two groups, those with g(r) < p (free or active cells), and 
those with g(r) > p (blocked or inactive cells). Denoting by q the density of inactive cells 
on the lattice we have q = 1 — p for < p < 1 and q = for p > 1. In this model the critical 
pressure is p c = 0.461. Periodic boundary conditions are used. We consider the evolution of 
the height of the i-th site in this model. Let us denote by hi(t) the height of the i-th generic 
site at time t. The set {hi,i = 1,...,N} where N = L/a, defines the interface between 
wet and dry cells. Given a site, chosen between N, say the site j, the height in the site i is 
increased by a with probability (i) 1 if j — i ± 1 and hi±\ > hi + 2a and hi < hi±2, (ii) 1/2 if 
j — i ± 1 and hi±\ > hi + 2a and hi = hi±2, (iii) 1 if j — i and hi < min(/ij_i, /ij+i) + 2a and 
Fi{hi + a) = 1. Otherwise no growth happens. Fi{hi + a) = Q(p — giihi + a)) is called the 



activity function |T3[ and O(x) is the unit step function defined as O(x) = 1 for x > and 
equals to otherwise, p is the microscopic driving force and giihi + a) is the quenched noise 
just above the interface distributed in the interval [0, 1]. Notice that the activity function F 
is the competition between the driving force and the quenched noise, so F is also a "noise" . 
Provided that the system size is large and that the intrinsic fluctuations are not too large 
0, the evolution equation for the height in a site i, in a short lapse r, is 

dhi a _ 

af = -G. + *. (2) 

where r is the mean lapse between successive election of any site and Gi are the micro- 



scopic growing rules for the evolution of the height at this site due that a site j is chosen at 
time t. Here rji is a Gaussian "thermal" noise with zero mean and covariance 

( , m(t)Vj(f}) = ^rG i 6 ij 6(t-irj. (3) 

Notice that in the notation of 0] the transition rate from from a configuration H to another 
H is W(H, H ) = i J2k 'K^fc? h>k + a )Gk Hj^k ^(X > hj)- So the first moment is ^ J2h' i^k ~ 
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hk)W(H, H') = ^Gk- As a consequence of the fact that subsequent configurations differ 
only in the height at one site all the moments are diagonal and proportional to the first 
moment §|. 

For this model 1T3IJTT 



Gi(hi-i, h h h i+1 ) = W i+1 + + Fi(\hi\ +l)Wi, (4) 

where \hi\ = [— ] denotes the integer part of hi in units of a. This definition is meaningless 
if a is taken as one, as in the discrete model. We shall show below that in the continuous 
limit, it means that F is taken as a constant inside a pore. In Eq. (|j), 

W i±1 = \ [1 - G(Hj ±2 ) + O(tfp)] ~ 2) , ~ 

Wi = l- G(Ui - 2) , {0) 

where H£ = (h s — h r )/a and Ui = -[hi — min(/i i+1 , Notice that all the heights are in 



units of a in order to keep the arguments of the step function without units. For Wi±\ [13 
the 5 Kronecker function has been taken as 

5{x,y) = Q{x-y) + e{y-x)-l. (6) 

Using the fact that min(x, y) = \{{x + y) — (x — y) [Q(x — y) — Q(y — x)]} and with a more 
compact notation 

Ui = H^m + Hi-i + H i-l - • (7) 

The representation of the step function can be expanded as Q(x) = Y,T=o °k xh providing 
that x is smooth. Our focus is on properties of the surface on large length scales so we 
kept the expansion of the step function to first order in his argument. The best choice 
for the representation of the function is the shifted hyperbolic tangent ]9[], defined as 
Q(x) = {1 + tanh[C(x + b)]}/2, where b is the shift and C is a parameter that allows to 
recover the in the limit C —>■ oo. We choose b = 1/2. The reason of our choices is that it 
allows us to define the S function as Eq. (|J). The coefficients fulfill 

1 c c c 

°o = g t 1 + tanh (^)} and ci = — cosh _2 (— ) . (8) 

We shall show bellow the qualitative information that can be obtained from these coefficients 
assuming this regularization scheme. The next step is to regularize the height defining an 
interpolating function. This is done by expanding the height h i+ i = h(xi + xi) around 
Xi = ia. Retaining only the leading terms in the expansion, the adimensional difference of 
heights is 

HttL = (i-m) d x h\ Xi + \{e - m 2 ) d 2 x h\ Xi a + 0{a 2 ) , (9) 
where d x h = d^h/dxK 

Notice that in any discrete model there is in principle an infinite number of nonlinear- 
ities, but at long wavelengths the higher order derivatives can be neglected using scaling 
arguments, since one expect affine interfaces over a long range of scales, and then one is 
usually concerned with the form of the relevant terms. 
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Replacing Eq. (§) in Eqs. (§) and (|7|), using the expansion of the step function and 
retaining the leading terms to order 0(a), Eq. (§) can be written as 

dk{ ^ = ~ [W{ Xi + a) + W( Xl - a) + W( Xi ) F(x h \h( Xi )\ + 1)] + V (x h t) , (10) 

with 

W(x + a) + W(x -a) = (c - 2ci) + Ac\ (d x h) 2 + ac x [± + 4(c - 2c 1 )] d 2 x h , (11) 

W(x) = 1 - (c - 2ci) - 4c 2 (^/i) 2 + ifl Cl <9 2 /i . (12) 

Notice that the argument of F = Q(p — g(x iy \h(xi) \ + 1)) is not smooth, so his expansion is 
meaningless. For the activity function F we define an interpolation function 

F{x u h{xi)) = F(x t , \hi\) + —[F(xi, \h\ + 1) - F(x u \hi\)] + 0(5h 2 ) , (13) 

(X 

with < 5h < a that measures the departure of the height from the low pore. Then F is 
a smooth function taking continuous values in the interval [0,1]. With this definition we 
ensure that the characteristic size of the correlation between pores is of the order of the 
pore size. In real materials there always exist a typical size of the inhomogeneities in the 
disordered media which plays the role of the lattice constant a. 

The final step is a coarse-grained spatial average of the variables in order to obtain 
smooth continuous functions at a macroscopic level. In this way we obtain the stochastic 
continuous equation for this model, 

dh 

°-g = fi(F) + u(F) %h + X(F) {d x hf + V (x, t) , (14) 



where 



MF) = [(c -2 Cl )(l-F) + F]-, (15) 

T 

v(F) = Cl [\ (1 + F) + 4(c - 2ci)] - , (16) 

T 

A(F)=4c 2 (l-F)-. (17) 

T 



and F = F(x, h) as was defined in Eq. (|T3"D. Notice that n(F) is now the effective competition 
between the driving force and the quenched noise. 

Equation ([14]) shows that the nonlinearity arises naturally as a consequence of the mi- 
croscopic model. 

This explains the previous numerical results obtained by Amaral et al. Q|, that studied 
the effects of an effective coefficient A e // from a tilted interface showing that the nonlinear 
term must exist. Our result is also in agreement with those of Reka et al. || that obtained 
numerically a parabolic shape of the local velocity as function of the gradient for the DPD 
model near above the criticality for different reduced forces (p/p c — 1). Now, let us con- 
sider for a moment the restrictions we have imposed over the coefficients Cq and c\ in the 
expansion of the representation of the function. When the conditions (0) are satisfied, the 
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coefficient v is always positive. The coefficient A is greater or equal to zero independently 
of the representation of the step function. On the other hand /x can be positive or nega- 
tive depending of the value of F. In Fig. [I] we show F as function of C for A = v (upper 
curve) and v = (lower curve). We can distinguish three regions. In the region \x < and 
A > v > 0, the function F takes values close to zero. Notice that in the discrete model, 
as we approach to the critical value F is mostly zero because the interface gets pinned by 
long chains of inactive sites jpftlTJJ. In the continuous model near the criticality \i mostly 



becomes negative braking the advance of the interface and F is close to zero. Then regions 
of small gradients will be halted by the effective driving force \x. In these case the main 
responsible of the nonlinearities is the lateral contribution [see Eq. (|TTj)] . This explain why 
this contribution enhances the roughness at the criticality as was predicted by Braunstein 
et al. @. 

Faraway above the criticality the nonlinear term becomes less relevant. In the limit 
p — > (p — > 1) we recover the KPZ (EW) equation with thermal noise. However the fact 
that these limits are recovered is a specific characteristic of this particular model. 

In Figs. we show the temporal scaling behavior of the roughness w obtained from the 
numerical integration of Eq. flli]) . At the criticality a slope (3 = 0.67 ± 0.05 was obtained. 
Above the threshold we recover a crossover between the exponent (3 = 0.67 ± 0.05 and the 
(3 m ~ 1/3 as was obtained by Leschhorn || by means of the numerical integration of the 
QKPZ equation and by his automaton version. In Fig. we show the scaling behavior 
of the correlation function C^r, t) = ([hi +r (t) — hi{t)] 2 ) 1 / 2 . The exponent obtained was 
a = 0.641 ± 0.07 in agreement with the DPD models. The numerical integration was made 
in short lattices. The results in large systems and the details of the integration will be 
published elsewhere. Notice that even if the exponents from our equation are very similar 
that the one obtained from the QKPZ one, our equation is very different. The main difference 
is that the coefficient of the nonlinear term in our equation is strongly affected by the local 
characteristic of the substratum. Another great difference is that the global velocity of the 
interface is affected by the coefficient of the Laplacian term. The study of these differences 
requires a great computational effort that is not in the aim of the present work. 

How our results could be used in order to explain the role played by the disordered media 
in the experiments? In the experiments the advancement of the interface is determinated 
by the coupled effect of the random distribution of the capillary sizes, the surface tension 
and the local properties of the flow, so it is not surprising that all these effect give rise to 
a multiplicative noise in any evolution equation that intends to represent an experimental 
growth with disordered media. 

Summarizing, we derive the continuous equation from the microscopic one for the TL 
model. In our work the nonlinear term arises naturally as consequence of the microscopic 
dynamics. The numerical integration of our equation reproduces very accurately the scaling 
exponents of the roughness. However these results show that the DPD model is not described 
by a QKPZ equation with additive noise even if the exponents of the roughness are very 
similar. Our equation also allows to explain that the lateral growth contribution is the main 
responsible of the roughness near the criticality. Finally, we hope that this framework can 
be used in other growing models with quenched noise. 
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FIGURES 



FIG. 1. F as function of C for p,(F) = (dashed line) and X{F) = u(F) (solid line). In the 
first region (I) X(F) > v(F) > and n(F) < 0, in the second region (II) X(F) > v(F) > and 
(x{F) > and in the third region (III) < X(F) < v(F) and y,{F) > 

FIG. 2. log-log plot of the square roughness w 2 vs time for C = 1.3. In (a) p = 0.1, for this 
value of C the critical pressure is p c — 0.1. In circles we shows the results obtained from the 
numerical integration of Eq. (jlj|). The dashed line is used as a guide and as exponent 2/3 = 1.34. 
In (b) p = 0.3, the dashed line has slope 2(5 = 1.34 and the solid line has slope 2/3 m = 0.66. The 
numerical integration has been done with L = 1024 and over 30 independent samples. 

FIG. 3. log-log plot of C*2(r) as function of r for p = 0.1 and C = 1.3. The dashed line that is 
used as a guide has slope 2a = 1.28 
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